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Abstract 

Continuous Galerkin Petrov time discretization scheme is tested on some 
Hamiltonian systems including simple harmonic oscillator, Kepler’s prob¬ 
lem with different eccentricities and molecular dynamics problem. In par¬ 
ticular, we implement the fourth order Continuous Galerkin Petrov time 
discretization scheme and analyze numerically, the efficiency and conser¬ 
vation of Hamiltonian. A numerical comparison with some symplectic 
methods including Gauss implicit Runge-Kutta method and general lin¬ 
ear method of same order is given for these systems. It is shown that the 
above mentioned scheme, not only preserves Hamiltonian but also uses 
the least CPU time compared with upto-date and optimized methods. 
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1 Introduction 

Non-dissipative phenomena arising in the fields of classical mechanics, molecular 
dynamics, accelerator physics, chemistry and other sciences are modeled by 
Hamiltonian systems. Hamiltonian systems define equations of motion based 
on generalised co-ordinates q, = (<71, <72, • •' 1 <7 n) and generalised momenta pi = 
{p\,P2, ■ ■ ■ ,Pn ) and are given as, 

dpi <)I I dqi dll 

dt dqi ’ dt dpi ’ 

having n degrees of freedom. H : R 2n x R 2ra —> R is the total energy 
Hamiltonian system. A separable Hamiltonian has the structure 

H(p,q) = T(p) + V(q) 


(1) 
of the 
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in mechanics, T = |P T M 1 P represents the kinetic energy and V being the 
potential energy. The Hamiltonian system in partitioned form takes the form 



The first observation is that, for autonomous Hamiltonian systems, H is an 
invariant, thus by differentiating H(p,q) with respect to time we have, 



We can write y = (p, q), then |1]) can be written as, 

l/ = J-'VH, 


where ' represents the derivative with respect to time, V is a gradient operator 
and J is a skew symmetric matrix consisting of zero matrix 0 and nxn identity 
matrix /, 


Another property of Hamiltonian systems is that its flow is symplectic, i.e. for 
a linear transformation d' : R 2n i—»• R 2rl , the jacobian matrix ^'(y) satisfies 


*' T (y)W(y) = J. 


Conservation laws for Hamiltonian systems are generally lost while integrat¬ 
ing these system. It is generally desirable to preserve the underlying qualitative 
property of solutions of Hamiltonian systems. This is achieved by using symplec¬ 
tic integrators from the class of one step, multistep and general linear methods. 
A lot of attention has been paid on the construction and implementation of such 
integrators, for details see Q], EL 0 and 0. 

The continuous Galerkin Petrov time discretization scheme (cGP) was in¬ 
vestigated in |5] for the system of ordinary differential equations (ODEs). In [B], 
this scheme was studied for the heat equation. In particular, the cGP(2) scheme 
has found to be 4th order accurate in the discrete time point and is A-stable 
method. 

The objective of this paper is to provide analysis of cGP(2) scheme 15} [6], 7 :) ;8] 
on some Hamiltonian systems and comparing it with other symplectic methods 
of order four including Gauss implicit Runge-Kutta method represented as irk4 
[9] and a g-symplectic general linear method represented by glm4 of same order 
developed in pm and m, In section two a brief introduction about the methods 
is given. The tested problems of Hamiltonian systems along with numerical 
experiments of these methods on Hamiltonian systems are described in third 
section. Conclusion based on numerical comparison of third section is given in 
fourth section. 
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2 The Methods 


2.1 Continuous Galerkin-Petrov method (cGP) 

As a model problem we consider the ODE system given in ©: Find u : [0, T m \ — > 
W such that 

d t u{t ) = F{t, u(t)) for f£(0,T m ), 

u(0 ) = 0 K ’ 

The weak formulation of problem © reads: Find u £ X such that u( 0) = uo 
and 

(F(t,u(t)),v(t)) dt V v £ Y, (3) 

where X denotes the solution space and Y the test space. To describe the time 
discretization of problem © let us introduce the following notation. We denote 
by I = [0,T m ] the time interval with some positive final time T m . We start 
by decomposing the time interval I into N subintervals I n := (t n —i,t n ), where 
n e {1,..., N} and 


(d t u(t),v(t)) dt = 


0 = to < *1 < • ■ ■ < tjv-l < tjv = Tra¬ 
in our time discretization, we approximate the continuous solution u(t) of prob¬ 
lem © on each time interval I n by a polynomial function: 

k 

u(t) « Uh(t ) := ^2 U 3 n (j) n j{t) V t £ I n , (4) 

3=0 

where the ’’coefficients” U-f are elements of the Hilbert space W and the basis 
functions <p n j £ P^(/ n ) are linearly independent elements of the standard space 
of polynomials on the interval I n with a degree not larger than a given order k. 

For a given time interval J C 1 and a Banach space B 1 we introduce the 
linear space of H-valued time polynomials with degree of at most k as 

P fc (J, B) := J u : J B : u{t ) = ^ U j t j , Vt G J, U 3 £ B, Vj 

[ 3=0 

Now, the discrete solution space for the global approximation Uh '■ I —> W is the 
space Xfi C A" defined as 

Xh i u £ C(I, W) : u\ In eF k (I n ,W) V n = 1,..., N} 

and the discrete test space is the space Yfi C Y given by 

Yh := {w £ L 2 (I, W) : u| Jb £ V k -ii!n, W) V n = 1,..., N}. 

The symbol h denotes the discretization parameter which acts in the error esti¬ 
mates as the maximum time step size h. := maxi<„<jv h n , where h n := t n —t n -\ 
is the length of the n-th time interval 

Let us denote by Xfc 0 := X]f D A 0 the subspace of Xft with zero initial 
condition. Then, it is easy to see that the dimensions of the spaces X£ 0 and 
Y^ coincide such that it makes sense to consider the following discontinuous 
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Galerkin-Petrov discretization of order k for the weak problem J3]) : Find Uh £ 
uo + Xfc 0 such that 

[ (d t u h (t),v h (t)) dt = f {F(t,u h (t)),v h (t))dt \/v h £Y£. (5) 

Jo Jo 

We will denote this discretization as the ” exact cGP(k )-method!\ Since the 
discrete test space Yjf is discontinuous, problem © can be solved in a time 
marching process. Therefore, we choose test functions Vhit) = vip n ^{t) with an 
arbitrary v € W and a scalar function 'tp n .i : I — > R which is zero on I \ I n and 
a polynomial ip n ^ £ Pfc_i(/ n ) on the time interval I n = [t n -i,t n \. Then, we 
obtain for each i = 0,..., k — 1 


(■ d t u h (t ), v) il> n ,i{t)dt = [ ( F{t , u h (t)), v) ip nii (t)dt V v £ W. 

Jin 


( 6 ) 


By the definition of the weak time derivative we get for Uh represented by © 
the equation 

k 

[ (dtUh{t),v) ip n ,i(t)dt = [ ([/%, v) „ (j)' n dt \/V £ W. 

Jin Jin J = 0 

We define the basis functions <j> n j £ P*,(/„) of © via the reference transforma¬ 
tion cu„ :/—>■/„ where I := [—1,1] and 

hr. 


t = w n (t) := 


tn—l T tn 


-t £ I n 


V t £ /, n = 1,..., N. 


Let (pj £ P fc (i), j = 0,..., k, be suitable basis functions satisfying the conditions 


^(-i) = $o,j, hi 1 ) = S k,j, 


(7) 


where Skj denotes the usual Kronecker symbol. Then, we define the basis 
functions on the original time interval I n by 


with t := cu m 1 (f) = — (t — 


tn, tn—\ 


£ I. 


Similarly, we define the test basis functions ip n> i by suitable reference basis 
functions xpi £ Pk-i(d), i.e., 

1pn,i(t) ■= W) V t £ I n , i = 0, . . . , k ~ 1. 

By the property ©, the initial condition and the continuity (with respect to 
time) of the discrete solution Uh : / —> W is equivalent to the conditions: 


U't = u 0 


and 


U° = U^_, V n > 2. 


We transform the integrals in ([6]) to the reference interval I and obtain the 
following system of equations for the ’’coefficients” £ W, j = 1,..., k, in the 
ansatz © : 


( u l’ v ) H = yI( F i Wn ^’Yl U nh(t) j 


dt 


V v £ W 

( 8 ) 
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where i = 0,..., k — 1, 


a it j := / d { $j(i)ipi(i) dt, 


and the ’’coefficient” 17° G IT is known. We approximate the integral on the 
right hand side of © by the (k + l)-point Gaufl-Lobatto quadrature formula: 

( u>n(t),^2 u i4>j(t) ] ,v\tj>i(i)di « \ F () > v 

\ 3=0 / / M=0 \ y 3= 0 / 



where are the weights and t^ £ [— 1 , 1 ] are the integration points with to = ~ 1 
and tk = 1. Let us define the mapped Gaufi-Lobatto points t n , M £ I n and the 
coefficients 7 by 

ln,/i -— /^ 2 ,/i. '— (7/x)) 73,M •— 7j(^M)• 

Then, the system © is equivalent to the following system of equations for the 
k unknown ’’coefficients” U-f £ W, j = 1,k, 


E“« = y ( F Vuew. ( 0 ) 

3=0 /i=0 \ \ 3=0 / / 

with the k ’’equations” i = 0 ,..., k— 1 where C/° = L r „_ 1 for n > 1 and 17° = uq- 
Once we have solved this system we enter the next time interval and set the 
initial value of the new time interval I n+ i to t /° +1 := Uf. If the Gaufi-Lobatto 
formula would be exact for the right hand side of © this time marching process 
would solve the global time discretization © exactly. Since in general there is 
an integration error we call the time marching process corresponding to © 
simply the ” cGP (k)- method'. 

In principle, we have to solve a coupled system for the U-f £ W which could 
be very expensive. However, by a clever choice of the functions <pj and ifi it 
is possible to uncouple the system to a large extend. In the following, we will 
discuss this issue for the special methods cGP(l), cGP(2) and for the general 
method cGP(fc), k > 3. In all cases, we choose the basis functions (f)j £ P*,(/) as 
the Lagrange basis functions with respect to the Gaufi-Lobatto points i.e., 

= ^ 7 > I 4 £ { 0 , • • •, fc}. 

Then, the method © reduces to 
k , k 

( u L v ) h = ( F ( tn d’ U h) , v ) Vv£W, i = 0,...,k- 1, 

3=0 3=0 

and by the choice of the test basis functions tfi £ Pfc_i(7) we try to get suitable 
values for the coefficients 07 , and /3ij. In the following, we will use the following 
abbreviation and assumption: 

Ff l (Ui):=F(t nJ ,Ui)£H' V j = 0,..., fc, n = 1,..., N. (10) 
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2.1.1 The cGP(l) method 

We use the 2-point Gaufi-Lobatto formula (trapezoidal rule) with wq = w\ = 1 
and to = — 1, ti = 1- The only test function ipo is chosen as i/jq (t) = 1. Then, 
we obtain 

Q!o,o = — 1, ol op = 1, /?o,o = A>,i = 1- 

Using the notation U n ~ 1 := Uh{t n -i) = U„ and U n := Uh(t n ) = U*, we obtain 
the following equation for the ’’unknown” U n £ W : 

(U n ,v) H - {U n ~\v) H = ^ {(P’(t n _ 1 ,^"- 1 ) +F(t n ,U n ),v)} 

for all v £ W which is the well-known Crank-Nicolson method. In operator 
notation it can be written in the equivalent form: 

u n = U n ~ l + y M_1 {Fitn-uU 71 - 1 ) + F(t n , U n )} . 


2.1.2 The cGP(2) method 

We use the 3-point Gaufi-Lobatto formula (Simpson rule) with wo = u >2 = 1/3, 
wi = 4/3 and t 0 = —1, t\ = 0, t 2 = 1. For the test functions A £ Pi(/), we 
choose 

(t) = 1. 

Then, we get 


(«id) 


M/2 1 -l/2\ 

l-i 0 1 ) 


(Ad) 


A/4 0 -l/4\ 

A/3 4/3 1/3 / 


and the assumption m, the system to compute the ”unknowns” U/, t/ 2 £ W 
from the known U° = , reads: 


Un = \u 0 n + \ul + ^M- x {F»{U° n )-Ft(Ul)} ( 11 ) 

U 2 n = UV+^-M-^F^U^ + AF^ + F^U*)}. ( 12 ) 

Let us denote the value for 17/ computed from m and depending on [7 2 by 
U* = G\(JJn) where G\ : W — > W in general is a nonlinear operator. We 
substitute this in the equation m and get, for the unknown [7 2 £ W, the 
following fixed point equation : 


U 2 n = G 2 n {Ul) := 17° + -fM- 1 {F° n {U° n ) + 4F/MM)) + F 2 (M} 


The mapping G 2 : W —> W is a contraction if the time step size x n is sufficiently 
small. 
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2.2 Gauss implicit Runge-Kutta methods 

For the general autonomous first order differential equations 


y'(t) = f(y(t)), 


(13) 


where for system (HD, we choose y 
Kutta methods are defined as 


P 

q 


and f(y) 


-V q V(q) 

-VpT(g) 


Runge- 


Vn+l = Vn + h Y bif (Yi) 
i= 1 

and 

s 

Yj = y„ + h aijf(Yj) 

i—l 

where the coefficients a^-, bi and stage s determine the method. The Gauss 
methods have the highest possible order r = 2s and are symplectic and symmet¬ 
ric. We exclusively consider s = 2, fourth order method for a fair comparison. 

2.3 General linear methods 

General linear methods provide numerical solutions of initial value problems of 
the form m A general linear method is of the form, 

Y = h(A ® I)f(Y) + {U® I)y [n ~ 1] , 
y [n] = h(B ® I)f(Y) + (V ® I)y [n ~ 1] . 

where A ® I is the Kronecker product of the matrix A and the identity matrix 
I and h represents the step size. The s— component vector Y are the stages 
and f(Y) are the stage derivatives. The vector j/b 1-1 ! with i —components is an 
input at the beginning of a step and results in output approximation ybd. With 
a slight abuse of notation, we can write, 

Y = hAf{Y) + Ui n ~ l \ 
yM = hBf(Y) + Vy [n ~ 1] . 

The matrices A , U, V and B represent a particular general linear method and 
are generally displayed as, 


' A 

U ' 

B 

V 


A fourth order symmetric G-symplectic general linear method is constructed 
with four stages (s = 4) and three input values (r = 3). The coefficeints of the 
method are given in uni. 


3 Numerical Experiments 

We performed numerical comparisons of the continuous Galerkin Petrov time 
discretization scheme , general linear method and implicit Gauss R-K method 
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all having the same order four, for some Hamiltonian systems including simple 
harmonic oscillator, Kepler’s problem with different eccentricities and molecular 
dynamical problems. Throughout the comparison, continuous Galerkin Petrov 
scheme is denoted by acronym cGP(2), while general linear method and implicit 
Gauss R-K method are represented by the acronym glm4 and irk4 respectively. 
The emphasis in our comparison is on the accuracy of solution, including the 
phase information, enrgy conservation and CPU time using above discussed 
methods. For each method and problem, we used different stepsizes and several 
intervals of integration. Stepsizes were chosen as a compromise between having 
small truncation error and performing efficient integration on each step. The 
accuracy of the solution was measured by the L 2 norm of the absolute global 
error in the position and velocity coordinates and is denoted by E g (t). The 
relative error in Hamiltonian is defined as 


E e (t) = 


m - m 

E( 0) 


Growth of global error is measured for first two problems as their exact solution 
exists, while relative error in Hamiltonian E e (t) is calcultaed for all problems. 
We also measured computational effort using the CPU time. All the comparisons 
are done on the same machine and are optimized using MATLAB. 


Simple Harmonic Oscillator 

As an example of simple harmonic oscillator a mass spring system having kinetic 
energy p 2 /(2m), where p = mv is the momentum of the system and potential 
energy ^ kq 2 . Where q is distance from the equilibrium, m is the mass of the 
body which is attached to spring and k is constant of proportionality often called 
as spring constant. Here the Hamiltonian is the total energy of the system and 
has one degree of freedom 


H( q , v ) = \k q ‘ + ^-. 

The equations of motion from the Hamiltonian are 

, 9H(q,p) , dH(q,p) 

q = -^- = P, P = -^-= ~q- 

op oq 

We compared the problem using different stepsizes of h = 0.005,0.01,0.025, 
and 0.05. Figure [T| gives the log-log graph for time versus global error E g (t) and 
relative error in Hamiltonian E e (t) using stepsize h = 0.005 for the time interval 
[0, 1000]. We found almost the same behavior of error growth for position and 
Hamiltonian using the rest of stepsizes. In Figure [U the top plot gives the 
growth of global error and is approximately same for all tested methods, irk4 and 
cGP(2) having the least error while glm4 with slightly bigger error. In bottom 
plot of figure |T] the error in Hamiltonian is conserved by the methods. We also 
calculated the error growth according to Brouwer’s law jT2], our calculation 
shows that the exponent of time is 1 and 0.6 for E g (t) and E e (t) respectively, 
closed to its expected value. Table |T| gives the cost of integration for simple 
harmonic oscillator using all stepsizes. The table lists the stepsizes, maximum 
of global error, maximum of Hamiltonian error and CPU time. We observe from 





Hamiltonian Error Global error 


Method 


stepsize ( h ) Max. of Global Max. of Hamiltonian CPU Time (sec.) 

Error Error 


cGP(2) 

0.05 

8.67 

X 

10~ 

-6 

4.08 

X 

IQ- 14 

2.8 

cGP(2) 

0.025 

5.42 

X 

10- 

-7 

1.07 

X 

10 -13 

5.3 

cGP(2) 

0.01 

1.38 

X 

10 - 

-8 

1.12 

X 

10“ 13 

12.1 

cGP(2) 

0.005 

8.68 

X 

10“ 

10 

1.31 

X 

10“ 13 

23.9 

glm4 

0.05 

2.51 

X 

10 - 

-5 

4.67 

X 

10- 11 

3.2 

glm4 

0.025 

1.57 

X 

10 - 

-6 

7.57 

X 

10- 13 

11.6 

glm4 

0.01 

4.09 

X 

10 - 

-8 

3.1 

X 

IO" 15 

78.6 

glm4 

0.005 

3.34 

X 

10 - 

-9 

1.14 

X 

10“ 15 

395.3 

irk4 

0.05 

8.67 

X 

10 - 

-6 

6.43 

X 

10“ 15 

5.7 

irk4 

0.025 

5.41 

X 

10 - 

-7 

2.51 

X 

10“ 14 

11.5 

irk4 

0.01 

1.31 

X 

10 - 

-8 

1.97 

X 

IQ- 14 

45.0 

irk4 

0.005 

6.98 

X 

10“ 

11 

3.68 

X 

IQ- 14 

191.8 


Table 1: Maximum of global error, Hamiltonian error and CPU time for simple 
harmonic oscillator. 


the Tablc[l]that cGP(2) used the least CPU time and also having the least value 
for maximum of global error except for h = 0.005, where irk4 having the least 
end point global error, may be because of entering in a dip also depicted in 
Figure [T| The methods irk4 and glm4 are using eight and sixteen times more 
CPU time than cGP(2) giving similar accuracy for h = 0.005. 


-cGP(2) 

-glm4 

1 irk4 


0 

0 

<3 

0 

o__ 

o> 

0 

: I -cGP(2) 

r- -glm4 

! | - irk4 

Tl .11 

} 


10 3 10 2 10 1 10° io 1 io 2 io 3 

Time 


Figure 1: The growth of global error and relative error in Hamiltonian for Simple 
harmonic oscillator using stepsize h = 0.005. 


Kepler’s Problem 

Kepler’s problem is two body orbital problem in which the bodies are moving 
under their mutual gravitational forces. We can assume that one body is fixed at 
the origin and the second body is located in the plane with coordinates (q±, 52 )- 
The solution of this problem is used in many important applications which 
includes the determination of orbits for new asteroids and the measurement 
of orbits for the two primary bodies in a restricted three body problem. The 
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Hamiltonian of the system can be written in separable form as [!_ 


This can be written as H = T+V, where T = (Pi+p|)/2 and V = — 1 / \Jq\ + q 2 
are kinetic and potential energy of the system respectively. As like the previous 
problem, this system is also autonomous so the Hamiltonian H is a conserved 
quantity. 

The equations of motion are 


<h = Pi, 
v'i = q" = 


v'i = <h 


q '2 =P2 
9i 

\ql + ql)* 

02 


with the initial conditions 


(14) 


0 i( 0 ) = 1 — e, 92 ( 0 ) = 0 , q[( 0 )= 0 , 92 ( 0 ) = 

where e is eccentricity 0 < e < 1. The exact solution of the above equations 
(fl4l) is 

yi = cos (E) — e, yi = \J\ — e 2 sin (E), 

and 


y' x = — sin(£')(l — ecos(E)) 1 , y' 2 = ^/(l — e 2 ) cos(m)(l — ecos (E)) 1 , 


where the eccentric anomaly E satisfies Kepler’s equation t = E — esin (E). 
Since Kepler’s equation is implicit in E , the equation is usually solved using 
a non-linear equation solver, although useful analytical approximations can be 
found for smaller eccentricity. 

The integrations are performed for Kepler’s problem with different eccen¬ 
tricities e = 0, 0.5 and 0.9. The integration is done for 1000 periods for e = 0 
and 100 periods for e = 0.5 and 0.9. For each method, we measured E g {t) and 
E e (t) throughout the interval of integration. A variety of different stepsizes are 
used to analyze the behaviour of error growth. We used the stepsizes of h = 
h = h = i§|q, h = and h = for eccentricities e = 0, 0.5 and 0.9. A 
log-log plot of time against error is given for Kepler’s problem in Figures [2] and 
fusing eccentricities 0 and 0.9 respectively. Growth of errors in both quantities 
behave in the same manner as for e=0.5. It is seen that the global error growth 
is approximately linear for cGP(2), irk4 and glm4, i. e., growing as t 0 9 (see 
figures [H and [3]). The error in Hamiltonian remains conserved for cGP(2), irk4 
and glm4 for the intervals of integration. Our calculation shows that for E e (t) 
grows as t °- 6 , showing a good agreement to its expected value. The cGP(2) ex¬ 
hibits a smaller error even the problem becomes more eccentricitic (see Figures 
[2] and [3| ■ 

We also measured the cost of integration for Kepler’s problem using all 
stepsizes for all three eccentricities. Tables EJ |3] and E] lists the stepsizes, 
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Hamiltonian Error Global Error 


Method 


stepsize ( h ) Max. of Global 

Error 


cGP(2) 

2tt/400 

4.88 x 10 -6 

cGP(2) 

271-/800 

2.54 x 10“ 7 

cGP(2) 

271-/1600 

6.56 x 10“ 8 

cGP(2) 

2tt/3200 

1.38 x 10 -8 

cGP(2) 

2tt/6400 

4.54 x 10“ 9 

glm4 

2tt/400 

2.36 x 10“ 5 

glm4 

27r/800 

1.56 x 10 -6 

glm4 

2tt/1600 

1.66 x 10 -7 

glm4 

2tt/3200 

8.59 x 10 -8 

glm4 

2tt/6400 

1.05 x 10 -8 

irk4 

2tt/400 

1.04 x 10“ 5 

irk4 

2tt/800 

5.68 x 10" 7 

irk4 

2tt/1600 

2.98 x 10“ 7 

irk4 

2tt/3200 

6.58 x 10 -8 

irk4 

2tt/6400 

1.27 x 10“ 8 


Max. of Hamiltonian CPU Time (sec.) 
Error 


4.07 

X 

10“ 

-13 

93.4 

7.54 

X 

10- 

-12 

192.7 

1.28 

X 

10- 

-11 

378 

2.28 

X 

io- 

-12 

760.5 

1.49 

X 

io- 

-12 

1487 

2.35 

X 

io- 

-14 

3464 

3.06 

X 

io- 

-14 

12172 

9.39 

X 

io- 

-14 

47724 

5.32 

X 

io- 

-13 

180716 

9.98 

X 

io- 

-12 

752864 

4.72 

X 

io- 

-14 

1598 

2.79 

X 

io- 

-14 

6348 

3.28 

X 

io- 

-14 

23609 

5.17 

X 

io- 

-13 

93215 

7.18 

X 

10- 

-12 

39460 


Table 2: Maximum of global error, Hamiltonian error and CPU time for Kepler’s 
Problem with e = 0 for 10 3 periods. 


maximum of global error, maximum of Hamiltonian error and CPU time for 
e = 0, e = 0.5 and 0.9 respectively. We observe from the information depicted 
in tables, that cGP(2) used the least CPU time and also having the least value 
for maximum of global error for all the stepsizes. For e=0, using the least step- 
size i.e h = g|^Q, irk4 and glm4 used 506 and 26 times more CPU time than 
cGP(2). While for e=0.5 and 0.9, irk4 and glm4 used nearly 55 and 24 times 
more CPU time than cGP(2). 



Time 


Figure 2: The growth of global error and relative error in Hamiltonian for 
kepler’s problem with e = 0 using stepsize 27r/6400 for 10 3 periods. 
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Method 


stepsize ( h ) Max. of Global Max. of Hamiltonian CPU Time (sec.) 

Error Error 


cGP(2) 

2?r/400 

1.06 x 10" 4 

2.83 x 10“ 8 

10.2 

cGP(2) 

2tt/800 

6.63 x 10~ 6 

1.77 x 10“ 9 

19.9 

cGP(2) 

2tt/1600 

4.15 x 10“ 7 

1.11 x 10~ 10 

39.8 

cGP(2) 

2tt/3200 

2.54 x 10“ 8 

7.01 x 10- 12 

79.2 

cGP(2) 

2tt/6400 

2.9 x 10' 9 

1.02 x 10- 12 

157.5 

glm4 

2tt/400 

2.1 x 10 -4 

5.97 x 10" 8 

43.8 

glm4 

2tt/800 

1.31 x 10“ 5 

3.74 x 10“ 9 

182.3 

glm4 

2tt/1600 

8.2 x 10~ 7 

2.34 x 10 _1 ° 

688.2 

glm4 

2tt/3200 

3.31 x 10“ 8 

1.46 x 10” 11 

2366 

glm4 

2tt/6400 

2.16 x 10 -8 

1.38 x 10- 12 

8789 

irk4 

2tt/400 

8.84 x 10 -5 

1.89 x 10 -8 

22.1 

irk4 

2tt/800 

5.53 x 10" 6 

1.18 x 10“ 9 

68.5 

irk4 

2tt/1600 

3.43 x 10“ 7 

4.71 x 10~ n 

244 

irk4 

2tt/3200 

1.04 x 10 -8 

4.66 x 10- 12 

955 

irk4 

2tt/6400 

2.51 x 10“ 8 

3.13 x 10- 13 

3830 


Table 3: Maximum of global error, Hamiltonian error and CPU time for Kepler’s 
Problem with e = 0.5 for 10 2 periods. 


Method 


stepsize (h) Max. of Global Max. of Hamiltonian CPU Time (sec.) 


Error 


cGP(2) 

2?r/400 

4.87 

cGP(2) 

2tt/800 

1.68 

cGP(2) 

2tt/1600 

1.92 x lO^ 1 

cGP(2) 

2tt/3200 

1.3 x 10~ 2 

cGP(2) 

2tt/6400 

8.41 x 10“ 6 

glm4 

2tt/400 

111.8 

glm4 

2tt/800 

4.54 

glm4 

2tt/1600 

1.46 

glm4 

2tt/3200 

9.77 x 10~ 2 

glm4 

2tt/6400 

6.14 x 10“ 3 

irk4 

2tt/400 

4.54 

irk4 

2tt/800 

1.27 

irk4 

2tt/1600 

1.32 x lO” 1 

irk4 

2tt/3200 

9.01 x 10“ 3 

irk4 

2tt/6400 

5.73 x 10“ 4 


Error 


5.23 

X 

10" 3 

10.4 

2.43 

X 

10 -4 

20.2 

1.42 

X 

10~ 5 

41.1 

8.74 

X 

io- 7 

80.6 

5.43 

X 

10 -8 

162.1 

6.5 : 

X 

10 -4 

44.1 

2.23 

X 

10- 4 

172 

1.62 

X 

10" 5 

676 

1.04 

X 

10" 6 

2482 

6.53 

X 

10 -8 

8402 

1.73 

X 

10 -5 

22.4 

1.36 

X 

10 -5 

67.8 

1.38 

X 

10-® 

257.5 

9.45 

X 

10' 8 

1086 

6.03 

X 

10' 9 

3786 


Table 4: Maximum of global error, Hamiltonian error and CPU time for Kepler’s 
Problem with e = 0.9 for 10 2 periods. 
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Figure 3: The growth of global error and relative error in Hamiltonian for 
kepler’s problem with e = 0.9 using stepsize 27r/1600 for 10 2 periods. 


Molecular Dynamical Problem 

We consider the interaction of seven Argon atoms in two dimension, where one 
of the atom is centered by six atoms which are symmetrically arranged m 
The Hamiltonian for the molecular dynamics is written as PQ 

H («>p) = \ E ^ p ^ pi +E E Vi i ii® - ® ii 

i=1 ' 1 »=2 j =1 


where I4,(r) are potential functions. Here q* andpi are positions and generelized 
momenta for the atoms. And m* denotes the atomic mass of the zth atom. 

The equations of motion for the frozen Argon crystals are given as 


g?(t) = 


24ecr 6 

Hli 


7 

E 

3= 


(gj - Qi) 

Lllgj-gilli 


- 2cr 6 


(gj - g») 


I|gi-gi|l2 J 


i = 


where r = 2, m* = 66.34 x 10 _27 [kg], cr,,' = cr = 0.341 [nm] and £ = 

1.654028284 x 10 —21 [J]. Initial positions and initial velocities are taken in [nm] 
and [nm/sec] respectively PQ. 

In molecular dynamics, since much ineterst is emphasized on macroscopic 
quantities like Hamiltonian. So we also discussed only the energy conservation of 
atoms over an interval of length 2 x 10 5 [fsec] (If sec = 10~ 6 ). The experiments 
are done using the stepsizes of 0.5 fsec, 1 fsec, 2 fsec and 4 fsec. The graphical 
results are only shown for h = 0.5 x 10~ 6 [fsec] as the error growth using other 
stepsizes was approximately same. Figure [I] shows that the tested methods 
conserve the value of Hamiltonian H even though the conservation is of highly 
oscillatory, while the error in Hamiltonian for cGP(2) grows as t 0,7 . On the other 
hand, for irk4 and glm4 the exponent of time is 0.59 and 0.61 respectively. 
Table [5] gives the cost of integration for molecular dynamical problem using 
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Relative Error in Hamiltonian 


Method 

stepsize ( h ) 

Max. of Global 
Error 

CPU Time (sec.) 

cGP(2) 

4 

1.24 x 

10 -9 

658 

cGP(2) 

2 

6.15 

X 

10- 11 

1323 

cGP(2) 

1 

7.05 

X 

10“ 12 

2680 

cGP(2) 

0.5 

3.73 

X 

io- 13 

5821 

glm4 

4 

9.24 

X 

io- 11 

1772 

glm4 

2 

5.58 

X 

10“ 12 

4300 

glm4 

1 

3.53 

X 

10“ 13 

11367 

glm4 

0.5 

9.41 

X 

10- 14 

34591 

irk4 

4 

3.72 

X 

io- 11 

2060 

irk4 

2 

2.32 

X 

10“ 12 

4272 

irk4 

1 

1.59 

X 

10“ 13 

10020 

irk4 

0.5 

2.34 

X 

10- 14 

23596 


Table 5: Maximum of Hamiltonian error and CPU time for molecular dynamical 
problem for 2 x 10 5 [fsec]. 


all stepsizes. The table lists the stepsizes, maximum of Hamiltonian error and 
CPU time. It is observed from the Table [5] that cGP(2) used the least CPU time 
for all the stepsizes used but exhibiting slightly big maximum of Hamiltonian 
error. The methods irk4 and glm4 having almost the same error growth for the 
integrated interval. 



Figure 4: The growth of relative error in Hamiltonian using h = 0.5 x 10 6 [f sec] 
for molecular dynamical problem over an interval of 2 x 10 5 [fsec]. 


4 Summary 

We implemented and analyzed the cGP(2) for Hamiltonian systems such as 
harmonic oscillator, Kepler’s problem and molecular dynamical problem. The 
obtained results are also compared with symplectic methods irk4 and glm4. It 
is shown that the cGP(2) method conserves the hamiltonian as other tested 
symplectic methods do. Moreover, giving the efficiency approximately same 
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as other methods yield, cGP(2) uses marginally less CPU time than compared 
methods. 
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